package deformation.line;

import grid.GridPoint;
import handler.DeformationLine;
import handler.DeformationPoint;
import handler.LineHandler;

import java.util.List;

import deformation.AbstractLineDeformation;
import deformation.Deformation;
import deformation.PreCompute;

public class AffineLineDeformation extends AbstractLineDeformation<AffineLineDeformation.PreData>{

	//alpha=2 (closed formula in the appendix of the paper
	//more closed formulas?
	final double  alpha;
	public AffineLineDeformation(LineHandler lineHandler) {
		super(lineHandler);
		 alpha=2;
	}

	public AffineLineDeformation(LineHandler lineHandler, double alpha) {
		super(lineHandler);
		this.alpha = alpha;
	}

	class PreData implements PreCompute{
		//A is a 1x2 matrix
		double [] A1;
		double [] A2;
		//here the weights (in the script W_i) is a 2x2 matrix depending on the delta**
		double[] delta00;
		double[] delta01;
		double[] delta11;
	}
	

	private PreData computeDeltaWeights(GridPoint v, List<DeformationLine> baseLines){
		
		PreData precomputedData=new PreData();
		precomputedData.delta00=new double[baseLines.size()];
		precomputedData.delta01=new double[baseLines.size()];
		precomputedData.delta11=new double[baseLines.size()];
		
		
		//compute the delta00, delta01, delta11
		
		if (alpha==2){
			for (int i=0;i<baseLines.size();i++){
				DeformationPoint a=new DeformationPoint(baseLines.get(i).getStartPoint().x,baseLines.get(i).getStartPoint().y);
				DeformationPoint b=new DeformationPoint(baseLines.get(i).getEndPoint().x,baseLines.get(i).getEndPoint().y);
				//orthogonal of (a_i-v) times (a_i-b_i)^transposed
				double bigDelta=(v.y-a.y)*(a.x-b.x)+(a.x-v.x)*(a.y-b.y);
				double beta00=(a.x-v.x)*(a.x-v.x)+(a.y-v.y)*(a.y-v.y);
				double beta01=(a.x-v.x)*(v.x-b.x)+(a.y-v.y)*(v.y-b.y);
				double beta11=(v.x-b.x)*(v.x-b.x)+(v.y-b.y)*(v.y-b.y);
				double theta=Math.atan(((b.x-v.x)*(b.x-a.x)+(b.y-v.y)*(b.y-a.y))/((v.y-b.y)*(b.x-a.x)+(b.x-v.x)*(b.y-a.y)))
							-Math.atan(((a.x-v.x)*(a.x-b.x)+(a.y-v.y)*(a.y-b.y))/((v.y-a.y)*(a.x-b.x)+(a.x-v.x)*(a.y-b.y)));
				precomputedData.delta00[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))
						*(beta01/beta00 - beta11*theta/bigDelta)/(2*bigDelta*bigDelta);
				precomputedData.delta01[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))
						*(1 - beta01*theta/bigDelta)/(2*bigDelta*bigDelta);
				precomputedData.delta11[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))
						*(beta01/beta11 - beta00*theta/bigDelta)/(2*bigDelta*bigDelta);		
				//System.out.println(precomputedData.delta00[i]+ "..."+ precomputedData.delta01[i] +"..."+ precomputedData.delta11[i]);
			}
		}
		else if (alpha==1){
			for (int i=0;i<baseLines.size();i++){
				DeformationPoint a=new DeformationPoint(baseLines.get(i).getStartPoint().x,baseLines.get(i).getStartPoint().y);
				DeformationPoint b=new DeformationPoint(baseLines.get(i).getEndPoint().x,baseLines.get(i).getEndPoint().y);
				double r=(b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y); //(b-a)^T*(b-a)
				double s=2*((b.x-a.x)*(a.x-v.x)+(b.y-a.y)*(a.y-v.y)); //2*(b-a)^T*(a-v)
				double t=(a.x-v.x)*(a.x-v.x)+(a.y-v.y)*(a.y-v.y); //(a-v)^T*(a-v)
				double x=1;
				double first_int=(4*r*r+4*r*(s-t)+s*s)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-(2*r+s)*Math.log(x*(r*x+s)+t)
						+2*r*(x-1);		
				x=0;		
				double second_int=(4*r*r+4*r*(s-t)+s*s)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-(2*r+s)*Math.log(x*(r*x+s)+t)
						+2*r*(x-1);	
				precomputedData.delta00[i]=Math.sqrt(r)*
					(first_int-second_int)
					/(2*r*r);
				
				x=1;
				first_int=(2*r*s-4*r*t+s*s)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-(r+s)*Math.log(x*(r*x+s)+t)
						+2*r*x;		
				x=0;		
				second_int=(2*r*s-4*r*t+s*s)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-(r+s)*Math.log(x*(r*x+s)+t)
						+2*r*x;	
				precomputedData.delta01[i]=-Math.sqrt(r)*
					(first_int-second_int)
					/(2*r*r);
				
				x=1;
				first_int=(2*r*s-4*r*t+s*s)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-(r+s)*Math.log(x*(r*x+s)+t)
						+2*r*x;		
				x=0;		
				second_int=(2*s*s-4*r*t)*Math.atan((2*r*x+s)/Math.sqrt(4*r*t-s*s))
						/(Math.sqrt(4*r*t-s*s))
						-s*Math.log(x*(r*x+s)+t)
						+2*r*x;	
				precomputedData.delta11[i]=-Math.sqrt(r)*
					(first_int-second_int)
					/(2*r*r);
//				System.out.println(4*r*r+4*r*(s-t)+s*s);
//				System.out.println(Math.atan((2*r+s)/Math.sqrt(4*r*t-s*s))-Math.atan((s)/Math.sqrt(4*r*t-s*s)));
//				System.out.println(Math.log(t/(t+s+r)));
//				System.out.println((2*r+s));
//				precomputedData.delta00[i]=Math.sqrt(r)*(
//						((4*r*r+4*r*(s-t)+s*s)/(Math.sqrt(4*r*t-s*s)))*(Math.atan((2*r+s)/Math.sqrt(4*r*t-s*s))-Math.atan((s)/Math.sqrt(4*r*t-s*s)))
//						+(2*r+s)*(Math.log(t/(t+s+r)))
//						+2*r)		
//						/(2*r*r);
//				precomputedData.delta01[i]=-Math.sqrt(r)*(
//						((2*r*s-4*r*t+2*s*s)/(Math.sqrt(4*r*t-s*s)))*(Math.atan((2*r+s)/Math.sqrt(4*r*t-s*s))-Math.atan((s)/Math.sqrt(4*r*t-s*s)))
//						+(r+s)*(Math.log(t/(t+s+r)))
//						+2*r)		
//						/(2*r*r);
//				precomputedData.delta11[i]=Math.sqrt(r)*(
//						((2*s*s-4*r*t)/(Math.sqrt(4*r*t-s*s)))*(Math.atan((2*r+s)/Math.sqrt(4*r*t-s*s))-Math.atan((s)/Math.sqrt(4*r*t-s*s)))
//						+(s)*(Math.log(t/(t+s+r)))
//						+2*r)		
//						/(2*r*r);
						
			}
		}
		else if (alpha==0){
			for (int i=0;i<baseLines.size();i++){
				DeformationPoint a=new DeformationPoint(baseLines.get(i).getStartPoint().x,baseLines.get(i).getStartPoint().y);
				DeformationPoint b=new DeformationPoint(baseLines.get(i).getEndPoint().x,baseLines.get(i).getEndPoint().y);
				precomputedData.delta00[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))/3;
				precomputedData.delta01[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))/6;
				precomputedData.delta11[i]=Math.sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))/3;
			}
		}
		return precomputedData;
	}

	@Override
	public GridPoint calculateGridPoint(GridPoint v, PreData precomputed) {
		List<DeformationLine> movedLines = lineHandler.getMovedLines();
		
		//get qStar
		DeformationPoint qStar=getStar(movedLines, precomputed.delta00, precomputed.delta01, precomputed.delta11);
		
		//values of the calculated gridpoint
		double xValue=qStar.x;
		double yValue=qStar.y;
		
		for (int j=0;j<movedLines.size();j++){
			//compute c^ and d^
			DeformationPoint cHat=new DeformationPoint(movedLines.get(j).getStartPoint().x-qStar.x,movedLines.get(j).getStartPoint().y-qStar.y);
			DeformationPoint dHat=new DeformationPoint(movedLines.get(j).getEndPoint().x-qStar.x,movedLines.get(j).getEndPoint().y-qStar.y);
			xValue+=precomputed.A1[j]*(cHat.x)+precomputed.A2[j]*(dHat.x);
			yValue+=precomputed.A1[j]*(cHat.y)+precomputed.A2[j]*(dHat.y);
		}
		
		return (new GridPoint(xValue,yValue));
	}


	@Override
	public PreData preComputeGridPoint(GridPoint v) {
		List<DeformationLine> baseLines = lineHandler.getBaseLines();
		PreData precomputedData=computeDeltaWeights(v, baseLines);
		
		//compute pStar
		DeformationPoint pStar=getStar(baseLines, precomputedData.delta00, precomputedData.delta01, precomputedData.delta11);
		
		
		//compute the inverse of the sum term in the A_j 
		//i.e. inv(sum_{all i} (^a_i;^b_i)^t * W_i * (^a_i;^b_i))
		//we name the matrix M=(m11,m12;m21,m22)
		double m11=0;
		double m12=0;  //equals m21
		double m22=0;
		for (int i=0; i<baseLines.size();i++){
			DeformationPoint a=baseLines.get(i).getStartPoint();
			DeformationPoint b=baseLines.get(i).getEndPoint();
			// m22 and m11 are exchanged for the future inversion
			m22+=precomputedData.delta00[i]*(a.x-pStar.x)*(a.x-pStar.x)
					+2*precomputedData.delta01[i]*(b.x-pStar.x)*(a.x-pStar.x)
					+precomputedData.delta11[i]*(b.x-pStar.x)*(b.x-pStar.x);
			
			m12+=precomputedData.delta00[i]*(a.x-pStar.x)*(a.y-pStar.y)
					+precomputedData.delta01[i]*((a.x-pStar.x)*(b.y-pStar.y)+(b.x-pStar.x)*(a.y-pStar.y))
					+precomputedData.delta11[i]*(b.x-pStar.x)*(b.y-pStar.y);			
			
			m11+=precomputedData.delta00[i]*(a.y-pStar.y)*(a.y-pStar.y)
					+2*precomputedData.delta01[i]*(b.y-pStar.y)*(a.y-pStar.y)
					+precomputedData.delta11[i]*(b.y-pStar.y)*(b.y-pStar.y);
		}
		//invert 
		double det=m11*m22-m12*m12;
		m11/=det;
		m22/=det;
		m12=-m12/det;
		
		//compute the A_j
		precomputedData.A1 =new double[baseLines.size()];
		precomputedData.A2 =new double[baseLines.size()];
		for (int j=0;j<baseLines.size();j++){
			//ahat=a-pStar
			DeformationPoint aHat= new DeformationPoint(baseLines.get(j).getStartPoint().x-pStar.x,baseLines.get(j).getStartPoint().y-pStar.y);
			//bhat=b-pStar			
			DeformationPoint bHat= new DeformationPoint(baseLines.get(j).getEndPoint().x-pStar.x,baseLines.get(j).getEndPoint().y-pStar.y);
			precomputedData.A1[j]=((v.x-pStar.x)*m11+(v.y-pStar.y)*m12)
					*(precomputedData.delta00[j]*(aHat.x)+precomputedData.delta01[j]*bHat.x)
					+((v.x-pStar.x)*m12+(v.y-pStar.y)*m22)
					*(precomputedData.delta00[j]*(aHat.y)+precomputedData.delta01[j]*bHat.y);
			precomputedData.A2[j]=((v.x-pStar.x)*m11+(v.y-pStar.y)*m12)
					*(precomputedData.delta01[j]*(aHat.x)+precomputedData.delta11[j]*bHat.x)
					+((v.x-pStar.x)*m12+(v.y-pStar.y)*m22)
					*(precomputedData.delta01[j]*(aHat.y)+precomputedData.delta11[j]*bHat.y);
		}
		
		return precomputedData;
	}
	
	/*
	 * computes and returns qStar or pStar
	 */
	static final DeformationPoint getStar(List<DeformationLine> lines,double[] delta00, double[] delta01, double[] delta11){
		double denominator=0;
		double starX=0;
		double starY=0;
		for (int i=0;i<lines.size(); i++){
			starX+=lines.get(i).getStartPoint().x*(delta00[i]+delta01[i])+lines.get(i).getEndPoint().x*(delta01[i]+delta11[i]);
			starY+=lines.get(i).getStartPoint().y*(delta00[i]+delta01[i])+lines.get(i).getEndPoint().y*(delta01[i]+delta11[i]);
			denominator+=delta00[i]+2*delta01[i]+delta11[i];
		}
		return (new DeformationPoint(starX/denominator,starY/denominator));
	}
	

}
